Figure3

Author

Affiliation

Isaias Hernandez

Postdoc - INSERM U1138 Centre de Recherche des Cordeliers

To reproduce Panels a-g you need to download data from Hu Junyi, et al. Multi-omic profiling of clear cell renal cell carcinoma identifies metabolic reprogramming associated with disease progression. Nature Genetics (2024) available at Zenodo, the rds file Databases_signatures, and the Supplementary Table S12+S7 from the article.

To reproduce Panels h and i, you need to download the Processed spatial transcriptomic data from the GEO and then obtain the Starfysh (click here to see how to install) immune deconvolution values by following the python code in the “Multiple-slide immune deconvolution by starfysh” script. Gene signatures to use with Starfysh can be found here: ccRcc_TLS_gene_signatures_starfysh.

Panel a

Code

#Needed data from Zenodo:
#bulkRNA_matrix_TPM.csv
#LC_matrix.csv
#GC_matrix.csv

library(ggplot2)
library(ggsci)
plot_correlation<-function(vect_x,
                             vect_y,
                             method="spearman",
                             smooth_method=c("lm","loess")[1],
                             cols = NA,
                             density=F,
                             log_scale=F,
                             point_size=0.1,
                             ylims=NA,
                             labls=c("x","y")){
  interesction <- intersect(names(vect_x),names(vect_y))
  df <- data.frame(cbind(vect_x[interesction],vect_y[interesction]))
  test <- cor.test(x=df[,1],y=df[,2],method = method)
  rho <- test$estimate
  pval <- test$p.value
  if(!all(is.na(cols))){
    plot <- ggplot(df, aes(x=df[,1], y=df[,2],colour=cols))
  }else{
    plot <- ggplot(df, aes(x=df[,1], y=df[,2]))
  }
  if(density){
    plot <- plot + geom_density_2d_filled(alpha=0.4)
  }

  plot <- plot + geom_point(size=point_size)
  plot <- plot + geom_smooth(method=smooth_method)
  if(log_scale){
    if(!all(is.na(ylims))){
      plot<- plot + scale_y_continuous(trans = "log10",breaks = base_breaks(), limits=ylims,labels = prettyNum) + theme(panel.grid.minor = element_blank())
    }else{
      plot<- plot + scale_y_continuous(trans = "log10",breaks = base_breaks(),labels = prettyNum) + theme(panel.grid.minor = element_blank())
    }
  }
  # if(!all(is.na(ylims))){
  #   plot <- plot + ylim(ylims)
  # }
  plot <- plot+ labs(subtitle=paste("cor=",format(rho,digits =3), " p = ",format(pval,digits = 3),sep=" "),
         title=paste(labls[1]," - ", labls[2]),
         x=labls[1],
         y=labls[2])

  res <- list(rho,plot,pval)
  return(res)
}
coding_combat2_mat <- read.csv("./bulkRNA_matrix_TPM.csv",row.names = 1) ###matrix of bulk RNA-seq
dim(coding_combat2_mat)

################processing metabolism data#################

library(data.table)
library(NormalyzerDE)
library(dplyr)

LC_matrix <- fread("./LC_matrix.csv") 
dim(LC_matrix)
LC_matrix <- data.frame(LC_matrix[,16:165],row.names = LC_matrix$Metabolites)
dim(LC_matrix)

GC_matrix <- fread("./GC_matrix.csv")
GC_matrix <- data.frame(GC_matrix[,11:160],row.names = GC_matrix$`Metabolite name`)
dim(GC_matrix)


all(colnames(LC_matrix)==colnames(GC_matrix))

matrix_meta <- rbind(LC_matrix,GC_matrix)
matrix_meta <- log2(matrix_meta)

####CHeck if there is something related to GABA, Gamma-aminobutyric acid
head(matrix_meta[rownames(matrix_meta)=="Gamma-aminobutyric acid",1:3])
#                          C115_N    C24_N    D36_N
#Gamma-aminobutyric acid 3.453304 2.425521 2.453644

####YES!!!

############
tmp_dat_met<-as.data.frame(t(matrix_meta))
tmp_dat_met$Sample<-rownames(tmp_dat_met)

############
###Now RNA
library(GSVA)
matrix <- coding_combat2_mat
##Load manual signatures
##Load signatures
Manual_signatures<-read_excel("Supplementary_tables.xls", sheet = "Table S12") #Get from the article, Table S12
Manual_signatures<-as.data.frame(Manual_signatures[-1,-1])
Manual_signatures<-as.list(Manual_signatures)

geneset<-list("GABA_sinthesis"=Manual_signatures$GABA_synthesis,"GABA_receptor"=Manual_signatures$GABA_receptor)

res_h_rna <- gsva(as.matrix(matrix), geneset, method = "ssgsea", parallel.sz = 0, verbose = T)

tmp_dat_rna<-as.data.frame(t(res_h_rna))
tmp_dat_rna$Sample<-rownames(tmp_dat_rna)

tmp_dat_met<-tmp_dat_met[order(factor(tmp_dat_met$Sample, levels =tmp_dat_rna$Sample)),]

table(tmp_dat_rna$Sample==tmp_dat_met$Sample)

tmp1<-tmp_dat_rna[,"GABA_sinthesis"]
names(tmp1)<-tmp_dat_rna$Sample
tmp2<-tmp_dat_met[,"Gamma-aminobutyric acid"]
names(tmp2)<-tmp_dat_met$Sample
    
tmp_p<-plot_correlation(tmp1,tmp2,labls = c("GABA_sinthesis","Gamma-aminobutyric acid"))
Correlation_met_rna<-tmp_p[[2]] #Extract_plot
Correlation_met_rna
ggsave(Correlation_met_rna,filename = "Speraman_correlation_met_rna_gaba.pdf",width = 10,height = 8)

ggsave(Correlation_met_rna,filename = "Speraman_correlation_met_rna_gaba.png",width = 10,height = 8)

Panel b

Code

#Needed objects from Zenodo:
#normal_snRNA.rds #Single nuclei data for the Normal adjacent tissue (NATs)
#snRNA.rds #Single nuclei data for the tumoral tissue
#Manual_signatures from panel a

##Import normal data sn
library(dplyr)
library(ggplot2)
library(cowplot)
library(harmony)
library(Seurat)
library(ggplot2)
library(cowplot)
library(ggsci)
library(Rcpp)
library(RColorBrewer)
library(ggpointdensity)
library(symphony)
library(viridis)
options(stringsAsFactors = F)
snRNA_NATs <- readRDS("./normal_snRNA.rds")

celltype <- c(brewer.pal(11,"Set3")[c(1,3:8,10:11)],pal_aaas(alpha = 0.7)(10)[-2],brewer.pal(8,"Set2"),brewer.pal(11,"Paired"))#[13:1]

unique(snRNA_NATs$sample)
unique(snRNA_NATs$celltype2)

#2 samples
##Preprocess the same as the tumoral data
VlnPlot(snRNA_NATs, c("nUMI","nGene","percent.mito"), group.by = "sample",pt.size = 0)

snRNA_NATs <- SetIdent(snRNA_NATs, cells = NULL, 
                  value = snRNA_NATs@meta.data$celltype2)

snRNA_NATs <- NormalizeData(snRNA_NATs)

snRNA_NATs <- FindVariableFeatures(snRNA_NATs, nfeatures = 2000)
VariableFeaturePlot(snRNA_NATs)

snRNA_NATs <- ScaleData(snRNA_NATs, verbose = T)
snRNA_NATs <- RunPCA(snRNA_NATs, npcs = 50, verbose = FALSE)


PCAPlot(snRNA_NATs,cols=celltype)
ElbowPlot(snRNA_NATs, ndims = 50)
snRNA_NATs <- RunHarmony(snRNA_NATs, c("sample"),max.iter.harmony = 20)
snRNA_NATs <- FindNeighbors(snRNA_NATs, reduction = "harmony", dims = 1:50, do.plot = T)
snRNA_NATs <- FindClusters(snRNA_NATs, resolution = 1)

snRNA_NATs  <- RunUMAP(snRNA_NATs, reduction = "harmony", dims = 1:50, n.components = 2)

snRNA_NATs <- CellCycleScoring(snRNA_NATs, g2m.features = cc.genes$g2m.genes,
                        s.features = cc.genes$s.genes)

###Take out the cells marked as unknown
snRNA_NATs<-subset(snRNA_NATs,celltype2!="unknown")

##Proyect in UMAP expression of GABA + celltypes
###Calculate signature for GABA_synthesis etc...
geneset<-list("GABA_sinthesis"=Manual_signatures$GABA_synthesis,"GABA_receptor"=Manual_signatures$GABA_receptor)

##Hallmark 54 genes
snRNA_NATs <- AddModuleScore(object = snRNA_NATs, features = list(Manual_signatures$GABA_synthesis), name = "GABA_synthesis")
##GABA_receptor
snRNA_NATs <- AddModuleScore(object = snRNA_NATs, features = list(Manual_signatures$GABA_receptor), name = "GABA_receptor")

##DOtplot
###Panel b left
DotPlot(snRNA_NATs, features = c("GABA_synthesis1","GABA_receptor1"),cluster.idents = T,group.by = "celltype2") + RotatedAxis() + 
    scale_colour_gradient2(low = "dodgerblue2", mid = "grey90", high = "red",midpoint = 0)  +ylab("celltype2")+xlab("Signatures")+ guides(color = guide_colorbar(title = "Score"))

##save RDS
saveRDS(snRNA_NATs,"./snRNA_NATs.rds")


###For the tumoral:

#combined_matrix
snRNA <- readRDS("./snRNA.rds")

##Check object..
VlnPlot(snRNA, c("nUMI","nGene","percent.mito"), group.by = "sample",pt.size = 0)
#
snRNA$IMM_subtype<-paste0("IMM_",snRNA$class)
###Process data
snRNA <- SetIdent(snRNA, cells = NULL, 
                  value = snRNA@meta.data$celltype2)

snRNA <- NormalizeData(snRNA)

snRNA <- FindVariableFeatures(snRNA, nfeatures = 2000)
VariableFeaturePlot(snRNA)
snRNA <- ScaleData(snRNA, verbose = T)
snRNA <- RunPCA(snRNA, npcs = 50, verbose = FALSE)
celltype <- c(brewer.pal(11,"Set3")[c(1,3:8,10:11)],pal_aaas(alpha = 0.7)(10)[-2],brewer.pal(8,"Set2"),brewer.pal(11,"Paired"))#[13:1]

PCAPlot(snRNA)

ElbowPlot(snRNA, ndims = 50)
snRNA <- RunHarmony(snRNA, c("sample"),max.iter.harmony = 20)
snRNA <- FindNeighbors(snRNA, reduction = "harmony", dims = 1:50, do.plot = T)
snRNA <- FindClusters(snRNA, resolution = 1)

snRNA  <- RunUMAP(snRNA, reduction = "harmony", dims = 1:50, n.components = 2)

snRNA <- CellCycleScoring(snRNA, g2m.features = cc.genes$g2m.genes,
                        s.features = cc.genes$s.genes)

###Calculate signatures
DefaultAssay(snRNA) <- "SCT"
##Hallmark 54 genes
snRNA <- AddModuleScore(object = snRNA, features = list(Manual_signatures$GABA_synthesis), name = "GABA_synthesis")
##GABA_receptor
snRNA <- AddModuleScore(object = snRNA, features = list(Manual_signatures$GABA_receptor), name = "GABA_receptor")

snRNA$celltype3<-snRNA$celltype
snRNA$celltype3[snRNA$celltype2=="1_B"]<-"Bcells"
snRNA$celltype3[snRNA$celltype2=="2_Plasma"]<-"Plasma"

##save RDS
saveRDS(snRNA,"./snRNA.rds")

Panel c

Code

#Needed objects from panel b
#snRNA processed after running panel b
library(dplyr)
library(ggplot2)
library(cowplot)
library(harmony)
library(Seurat)
library(ggplot2)
library(cowplot)
library(ggsci)
library(Rcpp)
library(RColorBrewer)
library(ggpointdensity)
library(symphony)
library(viridis)
library(clustree)
library(cluster)
options(stringsAsFactors = F)
snRNA <- readRDS("./snRNA.rds")
Epi_snRNA <- snRNA[,snRNA$celltype2=="Malignant"]
#
#Epi_snRNA <- readRDS("./Epi_snRNA.rds")


##Process data
Epi_snRNA <- NormalizeData(Epi_snRNA)

Epi_snRNA <- FindVariableFeatures(Epi_snRNA, nfeatures = 2000)
VariableFeaturePlot(Epi_snRNA)
Epi_snRNA <- ScaleData(Epi_snRNA, verbose = T)

#Performed an usupervised clustering changing resolution parameters
dr_dim.use<-1:30
Epi_snRNA <- RunPCA(object = Epi_snRNA, npcs = 30, 
                       verbose = T)
  dims.use = dr_dim.use
##Get significant PCAs

ElbowPlot(Epi_snRNA, ndims = 30)

print("performing JackStraw to determine the best PCs to use for clustering and dimension reduction")
      #Determine the ‘dimensionality’ of the dataset: choosing the number of PCs to used for further analysis.
      #jackstraw
number_pcs<-30
Epi_snRNA = JackStraw(Epi_snRNA, dims = max(1:number_pcs))

Epi_snRNA = ScoreJackStraw(Epi_snRNA, dims = 1:number_pcs)
dims.use = Epi_snRNA@reductions$pca@jackstraw@overall.p.values
dims.use = dims.use[dims.use[, 2] < 0.05, 1] #Generates a list of the dimensions to use that are the most informative...
##Apparently all dimensions 1:30 are significant...
Epi_snRNA <- FindNeighbors(Epi_snRNA, dims = dims.use,reduction = "harmony",do.plot = T)

#Epi_snRNA <- FindNeighbors(Epi_snRNA, reduction = "harmony", dims = 1:50, do.plot = T)

#####Find optimal cluster resulution
print("Performing Analysis to determine the best resolution parameter for clustering")

pcs<-as.matrix(Epi_snRNA@reductions$pca@cell.embeddings)
d <- as.dist(1 - cor(t(pcs)))

res <- seq(0.1, 1.0, by = 0.1)
res<-c(0.025,0.05,res)
library(cluster)
clust_nb <- c()
 v<-c()
for (i in 1:length(res)) {
  tryCatch({
print(i)
  Epi_snRNA <- FindClusters(Epi_snRNA, verbose = FALSE, resolution = res[i])

  clust_nb <- c(clust_nb, length(levels(Epi_snRNA@meta.data[, c(paste0("RNA_snn_res.", res[i]))])))
 ##Get silhouette value
  sil <- silhouette(as.numeric(Epi_snRNA$seurat_clusters), d, Fun = mean)
gc()
v <- c(v, mean(sil[, 3]))
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}

Epi_snRNA@meta.data[] <- lapply(Epi_snRNA@meta.data, function(x) {

  if (is.factor(x)) as.numeric(as.character(x)) else x

})

names(clust_nb)<-res 
names(v)<-res 
library(clustree)
pclust = clustree(Epi_snRNA, prefix = "RNA_snn_res.", node_size = 10, node_alpha = 0.8)

ggsave(pclust , file=file.path("clustree.png"), device="png", width=9, height=12)        



###Plot the  silhouette results...     
plot(names(v),v)
plot(NULL, NULL, xlim = c(min(as.numeric(clust_nb)), max(as.numeric(clust_nb))), 
     ylim = c(-0.1, 0.5), xlab = "Number of Clusters using malignant Spots", 
     ylab = "")
rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], 
     col = "#EAE9E9", border = FALSE)
grid(col = "white", lty = 1, lwd = 1.5)
points(clust_nb,v, pch = 19, col = ggplot2::alpha("black", 
                                                              0.8), cex = 1.5)
lines(clust_nb, v,col = "#E51718", lwd = 2, lty = 4)
mtext("Average silhouette width", side = 2, line = 2, cex = 1.5, 
      col = "#224A8D", las = 3)
axis(side = 4, at = seq(0, 1, 0.2), labels = c("0.0", "0.2", 
                                               "0.4", "0.6", "0.8", "1.0"))
invisible(dev.copy2pdf(file = "Silhuette_clusters.pdf",  width = 5.5, height = 5))

saveRDS(d,"distance_pca_Epi_snRNA.rds")

###According to this the best resolution is: 0.1 which gives rise to 5 clusters

rm(d);gc()
res_final=0.1        

Epi_snRNA <- FindClusters(object = Epi_snRNA, resolution = res_final)


Epi_snRNA  <- RunUMAP(Epi_snRNA, reduction = "harmony", dims = dims.use)
##CHeck how it turns using pca
Epi_snRNA  <- RunUMAP(Epi_snRNA, reduction = "pca", dims = dims.use,reduction.name = "umap_pca")
##Visualize
p_tmp<-CellDimPlot(
  srt = Epi_snRNA, group.by = c("seurat_clusters"),
  reduction = "UMAP",palcolor = c("#fc772d", "#c6d5a8", "#b48f60", "#a75565", "skyblue")
  #,theme_use = "theme_blank"
) 
##Note: this UMAP does not have the name of the clusters yet, this was produced after (see panel d)
ggsave(filename = "UMAP_clusters_ccRCC.pdf",p_tmp,width = 12,height = 12)

#Save object
saveRDS(Epi_snRNA,"./Epi_snRNA.rds")

Panel d

Code

#Needed objects from panel c
#Epi_snRNA processed after running panel c
#Databases_signatures object
library(dplyr)
library(ggplot2)
library(cowplot)
library(harmony)
library(Seurat)
library(ggplot2)
library(cowplot)
library(ggsci)
library(Rcpp)
library(RColorBrewer)
library(ggpointdensity)
library(symphony)
library(viridis)
library(clustree)
library(cluster)
#
Epi_snRNA <- readRDS("./Epi_snRNA.rds")

###Calculate some needed signatures from Supplementary Table S7
Normal_kidney_corpuscle_cells_signatures<-read_excel("Supplementary_tables.xls", 
sheet = "Table S7",skip = 1)
Normal_kidney_corpuscle_cells_signatures<-as.list(Normal_kidney_corpuscle_cells_signatures)
##Normal_kidney_corpuscle_cells_signatures, 
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Normal_kidney_corpuscle_cells_signatures$TAL), name = "Normal_TAL")
##
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Normal_kidney_corpuscle_cells_signatures$PT), name = "Normal_PT")
##
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Normal_kidney_corpuscle_cells_signatures$DCT), name = "Normal_DCT")
###
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Normal_kidney_corpuscle_cells_signatures$CD), name = "Normal_CD")
##
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Normal_kidney_corpuscle_cells_signatures$CNT), name = "Normal_CNT")
###
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Normal_kidney_corpuscle_cells_signatures$DTL), name = "Normal_DTL")

###
library(SCP)

##########DEGs and pathway annotation
###Takes around 4 hours....be patient
Epi_snRNA <- RunDEtest(srt = Epi_snRNA, group_by = "seurat_clusters", fc.threshold = 1, only.pos = FALSE,assay = "RNA",BPPARAM = BiocParallel::MulticoreParam(workers = 5))
gc()
VolcanoPlot(srt = Epi_snRNA, group_by = "seurat_clusters")

DEGs_all <- Epi_snRNA@tools$DEtest_seurat_clusters$AllMarkers_wilcox
##Export as table

 DEGs_all<-DEGs_all %>%
   group_by(group1) %>%
   arrange(p_val_adj)%>%
 filter(avg_log2FC>=0) 
 
write.table(DEGs_all,file = "DEGs_all_clusters_ccrcc_sn.tab",append = F,sep = "\t",quote = F,row.names = F) 


DEGs <- DEGs_all[with(DEGs_all, avg_log2FC >= 0.2 & p_val_adj < 0.05), ]
table(DEGs$group1)

##Annotate pathways..
library(dplyr)


ht_DE_sn_clusters <- FeatureHeatmap(
  srt = Epi_snRNA, group.by = "seurat_clusters", features = DEGs$gene, feature_split = DEGs$group1,
  species = "Homo_sapiens", db = c("GO_BP", "KEGG", "Hallmark"), anno_terms = TRUE, 
  #feature_annotation_palcolor = list(c("gold", "steelblue"), c("forestgreen")),
  height = 5, width = 4,assay = "RNA"
)
print(ht_DE_sn_clusters$plot)  
##Heatmap annotated...

###GSEA
####Enrichment analysis(GSEA)
Epi_snRNA <- RunGSEA(
  srt = Epi_snRNA, group_by = "seurat_clusters", db = c("GO_BP", "KEGG", "Hallmark"), species = "Homo_sapiens",
  DE_threshold = "p_val_adj < 0.05",BPPARAM = BiocParallel::MulticoreParam(workers = 5)
)
##Plot
GSEAPlot(srt = Epi_snRNA, group_by = "seurat_clusters", plot_type = "comparison",db=c("GO_BP", "KEGG","Hallmark"),topTerm = 5,padjustCutoff = 0.05) 

##Cluster 3 is Immune related (probably closer to TLS or immune infiltrated zone...)

##Let's check other db
##Specific group..group_use = "
p_gsea_cluster_0<-GSEAPlot(Epi_snRNA, db = c("GO_BP", "KEGG","Hallmark"), group_by = "seurat_clusters", group_use = "0",topTerm = 8)  

##For 1
p_gsea_cluster_1<-GSEAPlot(Epi_snRNA, db = c("GO_BP", "KEGG","Hallmark"), group_by = "seurat_clusters", group_use = "1",topTerm = 8)  

##For 2
p_gsea_cluster_2<-GSEAPlot(Epi_snRNA, db = c("GO_BP", "KEGG","Hallmark"), group_by = "seurat_clusters", group_use = "2",topTerm = 8)  


##For 3
p_gsea_cluster_3<-GSEAPlot(Epi_snRNA, db = c("GO_BP", "KEGG","Hallmark"), group_by = "seurat_clusters", group_use = "3",topTerm = 8)  

p_gsea_cluster_3
##For 4
p_gsea_cluster_4<-GSEAPlot(Epi_snRNA, db = c("GO_BP", "KEGG","Hallmark"), group_by = "seurat_clusters", group_use = "4",topTerm = 8)  
p_gsea_cluster_4

#O: PT, GABA synthesis, GABA receptor, Omega oxidation, MHCII, Fatty acid oxidation, fatty acid degradation, amino acid metabolism
#1: CD, CNT, DCT, Angiogenesis, cell cycle, omega oxidation, cytokine response, TNF, Interferon
#2: CNT, little GABA synthesis, cell cycle, FAS_pentose, EMT, Hypoxia, pEMT, p53, Oxidative phosphorylation, EMT
#3: TAL, CD, DCT, Motzer_fao_AMPK, this is rather patient-specific
#4: GABA receptor, TAL, CD, DCT, Angiogenesis


####Final annotation of clusters...
Epi_snRNA$annoted_ccrcc_clusters<-ifelse(Epi_snRNA$seurat_clusters=="0","0:PT-like", ifelse(Epi_snRNA$seurat_clusters=="1","1:ProInflammatory",ifelse(Epi_snRNA$seurat_clusters=="2","2:p53-EMT", ifelse(Epi_snRNA$seurat_clusters=="3","3:Patient-specific", "4:PI3K")) ))

Epi_snRNA$annoted_ccrcc_clusters<-factor(Epi_snRNA$annoted_ccrcc_clusters,levels = c("0:PT-like","1:ProInflammatory","2:p53-EMT","3:Patient-specific","4:PI3K"))

#
colors_clusters_ccrcc_sn<-c("#fc772d", "#c6d5a8", "#b48f60", "#a75565", "skyblue")
names(colors_clusters_ccrcc_sn)<-levels(Epi_snRNA$annoted_ccrcc_clusters)

######Calculate some signatures to plot...
####For cluster 0:PT-like
#Databases_signatures
Databases_signatures<-readRDS("Databases_signatures.rds") #Available in the Data folder in this github
####For cluster 0:PT-like
#GOBP_MONOCARBOXYLIC_ACID_CATABOLIC_PROCESS
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$GOBP_MONOCARBOXYLIC_ACID_CATABOLIC_PROCESS), name = "Carboxyic_acid_catabolic")
#KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$KEGG_VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION), name = "VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION")
#KEGG_TRYPTOPHAN_METABOLISM
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$KEGG_TRYPTOPHAN_METABOLISM), name = "TRYPTOPHAN_METABOLISM")
#HALLMARK_FATTY_ACID_METABOLISM
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$HALLMARK_FATTY_ACID_METABOLISM), name = "FATTY_ACID_METABOLISM")


#####For cluster 1:ProInflammatory 
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$HALLMARK_TNFA_SIGNALING_VIA_NFKB), name = "TNFA_SIGNALING_VIA_NFKB")
#
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$HALLMARK_INTERFERON_GAMMA_RESPONSE), name = "IFN_gamma")
#
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$HALLMARK_IL6_JAK_STAT3_SIGNALING), name = "IL6_JAK_STAT3")

#####For cluster 2:p53-EMT
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$KEGG_OXIDATIVE_PHOSPHORYLATION), name = "OXIDATIVE_PHOSPHORYLATION")
#
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION), name = "Hallmark_EMT")
#
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$HALLMARK_P53_PATHWAY), name = "p53_pathway")

#####For cluster 3:Patient-specific
#

#####For cluster 4:PI3K
#
Epi_snRNA <- AddModuleScore(object = Epi_snRNA, features = list(Databases_signatures$KEGG_PHOSPHATIDYLINOSITOL_SIGNALING_SYSTEM), name = "PI3K_sginaling")
#

##Final DotPlot
signatures_toPlot_final<-c("GABA_synthesis1","GABA_receptor1","Normal_PT1","Carboxyic_acid_catabolic1","VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION1","TRYPTOPHAN_METABOLISM1","FATTY_ACID_METABOLISM1","TNFA_SIGNALING_VIA_NFKB1","IFN_gamma1","IL6_JAK_STAT31","Normal_CNT1","OXIDATIVE_PHOSPHORYLATION1","Hallmark_EMT1","p53_pathway1","Normal_TAL1","PI3K_sginaling1")


dotplot_clusters_sn_ccrcc<-GroupHeatmap(
    srt = Epi_snRNA,
    features = signatures_toPlot_final,
    #
    group.by = c("annoted_ccrcc_clusters"),
    #heatmap_palette = "YlOrRd",
    #cell_annotation = c("RECIST_R_nR"),
    #cell_annotation_palette = c("Dark2"),
    show_row_names = TRUE, row_names_side = "left",
    add_dot = TRUE, add_reticle = TRUE,flip = F, column_title_rot = 45, show_column_names = TRUE,row_title_rot = 45,
    column_names_rot = 45,height = 4.3,width = 1.5
    #cluster_rows = T,
    #cluster_columns = T
    # heatmap_palette = "magma"
    # group_palcolor  = list(c("#F98982","#09C1C6"))
)
print(dotplot_clusters_sn_ccrcc$plot)

Panel e and f

Code

##Needed objects from panel d
#Epi_snRNA processed after running panel d
##Needed objects from panel b
#snRNA_NATs processed after running panel b

library(dplyr)
library(ggplot2)
library(cowplot)
library(harmony)
library(Seurat)
library(ggplot2)
#
Epi_snRNA <- readRDS("./Epi_snRNA.rds")
snRNA_NATs <- readRDS("./normal_snRNA.rds")

##Merge the needed cell together
PT_normal_ccrcc_integrated<-merge(subset(Epi_snRNA,annoted_ccrcc_clusters!="3:Patient-specific"),subset(snRNA_NATs,celltype2=="PT"))

#Integrate using Harmony for faster computing...
PT_normal_ccrcc_integrated <- NormalizeData(PT_normal_ccrcc_integrated)

PT_normal_ccrcc_integrated <- FindVariableFeatures(PT_normal_ccrcc_integrated, nfeatures = 2000)
VariableFeaturePlot(PT_normal_ccrcc_integrated)

PT_normal_ccrcc_integrated <- ScaleData(PT_normal_ccrcc_integrated, verbose = T)
PT_normal_ccrcc_integrated <- RunPCA(PT_normal_ccrcc_integrated, npcs = 50, verbose = FALSE)

PCAPlot(PT_normal_ccrcc_integrated)

ElbowPlot(PT_normal_ccrcc_integrated, ndims = 50)    

PT_normal_ccrcc_integrated <- RunHarmony(PT_normal_ccrcc_integrated, c("sample"),max.iter.harmony = 20)
PT_normal_ccrcc_integrated <- FindNeighbors(PT_normal_ccrcc_integrated, reduction = "harmony", dims = 1:50, do.plot = T)
PT_normal_ccrcc_integrated <- FindClusters(PT_normal_ccrcc_integrated, resolution = 1)

PT_normal_ccrcc_integrated  <- RunUMAP(PT_normal_ccrcc_integrated, reduction = "harmony", dims = 1:50, n.components = 2)

##Visualize
PT_normal_ccrcc_integrated$annoted_ccrcc_clusters_PT<-ifelse(is.na(PT_normal_ccrcc_integrated$annoted_ccrcc_clusters),"normal-PT",PT_normal_ccrcc_integrated$annoted_ccrcc_clusters)

saveRDS(PT_normal_ccrcc_integrated,file = "PT_normal_ccrcc_integrated.RDS")


#Remove previous objects
rm(Epi_snRNA,snRNA_NATs)

######Perform trayectory analysis using SlingShot from SCP
library(SCP)

#Trajectory inference
##################
#RunSlingshot
PT_normal_ccrcc_integrated <- RunSlingshot(srt = PT_normal_ccrcc_integrated, group.by = "annoted_ccrcc_clusters_PT", reduction = "UMAP",start = "normal-PT")

FeatureDimPlot(PT_normal_ccrcc_integrated, features = paste0("Lineage", 1:2), reduction = "UMAP", theme_use = "theme_blank")


###THis is the panel e
CellDimPlot(PT_normal_ccrcc_integrated, group.by = "annoted_ccrcc_clusters_PT", reduction = "UMAP", lineages = paste0("Lineage", 1:2), lineages_span = 0.1)

##Dynamic features.... for panel f

PT_normal_ccrcc_integrated <- RunDynamicFeatures(srt = PT_normal_ccrcc_integrated, lineages = c("Lineage1", "Lineage2"), n_candidates = 200,BPPARAM = BiocParallel::MulticoreParam(workers = 15)
                                                   #,features = genes_to_keep2 #set to genes that are characteristic of each population...along with metaprogrammes (SMP)
                                                     )


DynamicPlot(
  srt = PT_normal_ccrcc_integrated, lineages = c("Lineage1", "Lineage2"), group.by = "annoted_ccrcc_clusters_PT",
  features = c("CA9","PAX8","GABA_sinthesis1"),
  compare_lineages = TRUE, compare_features = FALSE
)

Panel g

Code

#You need the spatial transcriptomics and metabolomics from Zenodo:
#This is an example for sample R29T
#Manual_signatures from Supplementary Table S12

spatial_correlation2<-function(spat_obj,var,threshold=0.5,cortest="pearson",colors_pal=c("firebrick1","gray72","palegreen","dodgerblue")){
    values <- NULL
    for(v in var){
      if(v %in% rownames(spat_obj)){
        spat_obj <- AddMetaData(spat_obj,
                                ifelse(spat_obj@assays$SCT@data[v,] > quantile(spat_obj@assays$SCT@data[v,],threshold),"high","low"),
                                col.name = paste0(v,"_cat"))
        values[[v]] <- spat_obj@assays$SCT@data[v,]
        
      }else{
        spat_obj <- AddMetaData(spat_obj,
                                ifelse(spat_obj@meta.data[,v] > quantile(spat_obj@meta.data[,v],threshold),"high","low"),
                                col.name = paste0(v,"_cat"))
        values[[v]] <- spat_obj@meta.data[,v]
      }
      
    }
    var_cat <- paste0(var,sep="_cat")
    spat_obj  <- AddMetaData(spat_obj,
                             factor(apply(spat_obj@meta.data[,var_cat],1,function(x) paste0(x,collapse = "_")),levels=rev(c("low_low",
                                                                                                                            "low_high",
                                                                                                                            "high_low",
                                                                                                                            "high_high"))),
                             col.name = paste0(var,collapse = "_"))
    get_cols <- which(table(spat_obj[[paste0(var,collapse = "_")]][,1])!=0)
    get_cols <- brewer.pal(4,"Spectral")[get_cols]
    # p1 <- SpatialDimPlot(spat_obj,
    #                      paste0(var,collapse = "_"),
    #                      cols = get_cols)
    p1 <- SpatialPlot(spat_obj,
                      paste0(var,collapse = "_"),
                      image.alpha = 0,
                      cols = get_cols)
    
    corres <- cor.test(values[[var[1]]],
                       values[[var[2]]],method = cortest,na.action="na.omit")
    
    
    table_cont <-spat_obj@meta.data[,var_cat]
    order_factor <- rev(c("high_high","high_low","low_high","low_low") )
    
    p2 <- plot_contingency_one_stack(data = table_cont,
                                     variables = colnames(table_cont),
                                     create_labels = F,
                                     order_factor = order_factor,
                                     palette = rev(brewer.pal(4,"Spectral")))
    
    p_cor <- ifelse(round(corres$p.value,digits = 3) == 0,"p < 10e-16",paste0("p = ",round(corres$p.value,digits = 3)))
    p_chisq <- ifelse(round(p2[[2]]$p.value,digits = 3) == 0,"p < 10e-16",paste0("p = ",round(p2[[2]]$p.value,digits = 3)))
    
    #p1 <- p1 + ggtitle(paste0(var,collapse = " "),
    #                   paste0("Correlation = ",round(corres$estimate,digits = 3)," ", p_cor))
    p1 <- p1 + ggtitle(paste0(var,collapse = " "),
                       paste0("Chisq.test: ",p_chisq, "\n Rho: ", round(corres$estimate,digits = 3)))+ scale_fill_manual(values=colors_pal) + DarkTheme() + hide_axis
    #####Change title to red and bold if significant chisq.test (<0.01)
    if (p2[[2]]$p.value<=0.01&corres$estimate>0) {
      p1 <- p1 +theme(plot.subtitle = element_text(color="red", face = "bold"))
    }
    
    #compute correlation plots
    v1 <- setNames(values[[1]],colnames(spat_obj))
    v2 <- setNames(values[[2]],colnames(spat_obj))
    
    corr <- plot_correlation(vect_x =v1,vect_y = v2,labls = var,method = cortest)
    p3 <- corr[[2]]+theme_classic()
    return(list("spatial"=p1,
                "cont_table"=p2,
                "spatial_cor"=p3))
  }

library(Cardinal)
library(ggpubr)
library(ggsci)
library(ggrepel)
########number II
#R29_T, high expression
R29_T_neg <- readImzML(name="R29-T-neg", folder = "spatial_metabolomics/") #Path were you have the downloaded the data

###This is for number (II) from panel g
image(R29_T_neg,mz=102.05623, normalize.image="none",colorscale=c(rep("black",2),col.map("jet",n=99)),contrast.enhance="histogram",smooth.image = "gaussian")

############Now the remaining numbers using spatial transcriptomics

visium_directory<-"Path_to_visium" #Path to download data from sample R29_T
seurat_obj <- Seurat::Load10X_Spatial(paste0(visium_directory,"/outs"))

#Low-quality / dying cells often exhibit extensive mitochondrial contamination
    seurat_obj[["percent.mito"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT.")
    #VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
    # Collect all genes coded on the mitochondrial genome
    mt.genes <- grep(pattern = "^MT-", x = rownames(seurat_obj), value = TRUE)
    #remove mt genes
    genes_to_keep <- setdiff(names(which(Matrix::rowSums(seurat_obj@assays$Spatial$counts )>min.genes)),mt.genes)#Keep genes with count superior to 5
    #filter cells passing threashold
    seurat_obj_before<-seurat_obj
    percent.mito_manual=30
    seurat_obj <- subset(seurat_obj,features =genes_to_keep, subset = nFeature_Spatial > 300 & percent.mito < percent.mito_manual)##Spatial spots featuring more than 30% of mitochondrial genes and less than 300 genes were filtered out, as they identified necrotic or damaged tissue areas which were validated by a pathologist (VV). Genes with counts in less than 5 spatial spots were discarded.
    cat("Spots removed: ", ncol(seurat_obj_before) - ncol(seurat_obj), "\n")
    cat("Genes kept: ", length(genes_to_keep),"from",nrow(seurat_obj), "\n") 
 

##Normalize expression data in Seurat..
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
      #identify highly variable genes
seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 2000)
      #Scale data prior to dimension reduction  ( stored in pbmc[["RNA"]]@scale.data )
all.genes <- rownames(seurat_obj)
vars.to.regress="percent.mito"
seurat_obj <- ScaleData(seurat_obj, features = all.genes,vars.to.regress = vars.to.regress)


###Calculate GABA synthesis and GABA receptor score
##Gaba synthesis as mean
intersect_tmp<-intersect(Manual_signatures$GABA_synthesis,rownames(seurat_obj))
print(length(intersect_tmp))

seurat_obj <- AddMetaData(seurat_obj,apply(as.matrix(seurat_obj[["Spatial"]]$data[intersect_tmp,]),2,mean),col.name = "GABA_synthesis_mean")

##Gaba receptor as mean
intersect_tmp<-intersect(Manual_signatures$GABA_receptor,rownames(seurat_obj))
print(length(intersect_tmp))

seurat_obj <- AddMetaData(seurat_obj,apply(as.matrix(seurat_obj[["Spatial"]]$data[intersect_tmp,]),2,mean),col.name = "GABA_receptor_mean")


####Import signatures from Supplementary Table S8 to ccrcc_sn_clusters_signatures

##Calculate PT signature from sn using the signature list...
intersect_tmp<-intersect(ccrcc_sn_clusters_signatures$`0:PT-like`,rownames(seurat_obj))
print(length(intersect_tmp))

seurat_obj <- AddMetaData(seurat_obj,apply(as.matrix(seurat_obj[["Spatial"]]$data[intersect_tmp,]),2,mean),col.name = "ccRCC_PT_like")


###################Plot for number (i), H&E image

he <- SpatialPlot(seurat_obj,repel = F,label = F,image.alpha=1,alpha = c(0,0)) + NoLegend()

###################Plot for number (iii), GABA_synthesis_mean spatially

#Use Spata
library(SPATA2)
###Extra for SPATA2
tmp_spata<-SPATA2::asSPATA2(
    object = seurat_obj,
    sample_name = tmp_name,
    image_name = "slice1",
    spatial_method = "Visium"
)

pl3<-SPATA2::plotSurface(object = tmp_spata,
            color_by = c("GABA_synthesis_mean"), display_image=T,#alpha_by = "SCGB3A1",
            smooth = F, normalize =F,#smooth_span = 0.5,
            pt_size = 1.2) + 
    ggplot2::labs(title = "GABA_synthesis_mean")
  



###################Plot for number (iv), ccRCC_PT_like spatially
pl4<-SPATA2::plotSurface(object = tmp_spata,
            color_by = c("ccRCC_PT_like"), display_image=T,#alpha_by = "SCGB3A1",
            smooth = F, normalize =F,#smooth_span = 0.5,
            pt_size = 1.2) + 
    ggplot2::labs(title = "ccRCC_PT_like")

###############Plot for number (v), the spatial correlation

tmp<-spatial_correlation2(spat_obj = seurat_obj,var = c("ccRCC_PT_like","GABA_synthesis_mean"),cortest  = "spearman")
tmp$spatial 

Panel h

Code

###To reproduce this panel you need:
#1: Download all the slides visium data from GEO
#2: Download the file "Signatures_Starfysh.csv" found in the "Data" folder
#3: Obtaine starfysh deconvolution values for each slide by Runnning the python script: "Multiple-slide immune deconvolution by starfysh"
#Manual_signatures

library(Seurat)


process_Visium_seurat<-function(seurat_obj_path,percent.mito_manual=30,min.genes=5,nFeature=300,normalization_method="SCTransform",number_pcs=50,use.jack=T,resolution_clus=1,Calculate_MCPcounter=T,generate_plots=T,plots_path=NULL,marker.test="roc",vars.to.regress="percent.mito",only.pos = TRUE, run_DE=T,...){
    raw_data_directory <-seurat_obj_path #This must be the "outs" from SpaceRanger
    #Load object
    seurat_obj <- Seurat::Load10X_Spatial(raw_data_directory)
    #Low-quality / dying cells often exhibit extensive mitochondrial contamination
    seurat_obj[["percent.mito"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT.")
    
    #VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
    
    # Collect all genes coded on the mitochondrial genome
    mt.genes <- grep(pattern = "^MT-", x = rownames(seurat_obj), value = TRUE)
    #remove mt genes
    genes_to_keep <- setdiff(names(which(Matrix::rowSums(seurat_obj@assays$Spatial$counts )>min.genes)),mt.genes)#Keep genes with count superior to 5
    #filter cells passing threashold
    seurat_obj_before<-seurat_obj
    seurat_obj <- subset(seurat_obj,features =genes_to_keep, subset = nFeature_Spatial > nFeature & percent.mito < percent.mito_manual)##Spatial spots featuring more than 30% of mitochondrial genes and less than 300 genes were filtered out, as they identified necrotic or damaged tissue areas which were validated by a pathologist (VV). Genes with counts in less than 5 spatial spots were discarded.
    cat("Spots removed: ", ncol(seurat_obj_before) - ncol(seurat_obj), "\n")
    cat("Genes kept: ", length(genes_to_keep),"from",nrow(seurat_obj), "\n") 
    
    cat("Performing: ", normalization_method," Normalization","\n") 
    ###Data normalization
    if(normalization_method=="SCTransform")
    {
      if (use.jack) {
        cat("You selected SCTransform and JackStraw, the following normalization is only to make JackStraw work\n it is not going to be used for other downstream analyses","\n") 
        #Do this only in case we need the JackStraw  since it only supports non-SCT assays
        seurat_obj_tmp_jack<-seurat_obj
        seurat_obj_tmp_jack <- NormalizeData(seurat_obj_tmp_jack, normalization.method = "LogNormalize", scale.factor = 10000)
        #identify highly variable genes
        seurat_obj_tmp_jack <- FindVariableFeatures(seurat_obj_tmp_jack, selection.method = "vst", nfeatures = 2000)
        #Scale data prior to dimension reduction  ( stored in pbmc[["RNA"]]@scale.data )
        all.genes <- rownames(seurat_obj_tmp_jack)
        seurat_obj_tmp_jack <- ScaleData(seurat_obj_tmp_jack, features = all.genes,vars.to.regress = vars.to.regress)
        
      }
      seurat_obj <- SCTransform(seurat_obj, assay = "Spatial", verbose = T)
    }else{
      #global scaling normalization
      seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
      #identify highly variable genes
      seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 2000)
      #Scale data prior to dimension reduction  ( stored in pbmc[["RNA"]]@scale.data )
      all.genes <- rownames(seurat_obj)
      seurat_obj <- ScaleData(seurat_obj, features = all.genes,vars.to.regress = vars.to.regress)
    }
    #
    #reduction, calculate PCA, you can manually change the number of pcs to use...
    seurat_obj <- RunPCA(seurat_obj, features = VariableFeatures(object = seurat_obj),npcs = number_pcs)
    
    ###Perform further analysis using the best dimensions or manually set dimensions
    if (use.jack==TRUE&normalization_method=="SCTransform") {
            print("performing JackStraw to determine the best PCs to use for clustering and dimension reduction")
      #Determine the ‘dimensionality’ of the dataset: choosing the number of PCs to used for further analysis.
      #jackstraw
      seurat_obj_tmp_jack <- RunPCA(seurat_obj_tmp_jack, features = VariableFeatures(object = seurat_obj_tmp_jack),npcs = number_pcs)
      seurat_obj_tmp_jack = JackStraw(seurat_obj_tmp_jack, dims = max(1:number_pcs))
      seurat_obj_tmp_jack = ScoreJackStraw(seurat_obj_tmp_jack, dims = 1:number_pcs)
      dims.use = seurat_obj_tmp_jack@reductions$pca@jackstraw@overall.p.values
      dims.use = dims.use[dims.use[, 2] < 0.05, 1] #Generates a list of the dimensions to use that are the most informative...
      seurat_obj <- FindNeighbors(seurat_obj, dims = dims.use)
      seurat_obj <- FindClusters(seurat_obj, resolution = resolution_clus)
      ##Calculate also phenograph clusters..
      #Other clustering method:  https://github.com/sararselitsky/FastPG
      cluster = FastPG::fastCluster(seurat_obj@reductions$pca@cell.embeddings[,dims.use], num_threads = 35)
      seurat_obj$cluster_phenograph = cluster$communities
      #Run UMAP
      seurat_obj <- RunUMAP(seurat_obj, dims = dims.use)
      ###Run TSNE
      seurat_obj <- RunTSNE(seurat_obj, reduction = "pca", dims = dims.use,do.fast = T, k.seed = 10, check_duplicates = FALSE,perplexity = 30)
      
    }
    if (use.jack==TRUE&normalization_method!="SCTransform") {
      
      print("performing JackStraw to determine the best PCs to use for clustering and dimension reduction")
      #Determine the ‘dimensionality’ of the dataset: choosing the number of PCs to used for further analysis.
      #jackstraw
      seurat_obj = JackStraw(seurat_obj, dims = max(1:number_pcs))
      seurat_obj = ScoreJackStraw(seurat_obj, dims = 1:number_pcs)
      dims.use = seurat_obj@reductions$pca@jackstraw@overall.p.values
      dims.use = dims.use[dims.use[, 2] < 0.05, 1] #Generates a list of the dimensions to use that are the most informative...
      seurat_obj <- FindNeighbors(seurat_obj, dims = dims.use)
      seurat_obj <- FindClusters(seurat_obj, resolution = resolution_clus)
      ##Calculate also phenograph clusters..
      #Other clustering method:  https://github.com/sararselitsky/FastPG
      cluster = FastPG::fastCluster(seurat_obj@reductions$pca@cell.embeddings[,dims.use], num_threads = 35)
      seurat_obj$cluster_phenograph = cluster$communities
      #Run UMAP
      seurat_obj <- RunUMAP(seurat_obj, dims = dims.use)
      ###Run TSNE
      seurat_obj <- RunTSNE(seurat_obj, reduction = "pca", dims = dims.use,do.fast = T, k.seed = 10, check_duplicates = FALSE,perplexity = 30)
      
    }
    if (use.jack==FALSE)
      {
      seurat_obj <- FindNeighbors(seurat_obj, dims = 1:number_pcs)
      seurat_obj <- FindClusters(seurat_obj, resolution = resolution_clus)
      ##Calculate also phenograph clusters..
      #Other clustering method:  https://github.com/sararselitsky/FastPG
      cluster = FastPG::fastCluster(seurat_obj@reductions$pca@cell.embeddings[,1:number_pcs], num_threads = 35)
      seurat_obj$cluster_phenograph = cluster$communities
      #Run UMAP
      seurat_obj <- RunUMAP(seurat_obj, dims = 1:number_pcs)
      ###Run TSNE
      seurat_obj <- RunTSNE(seurat_obj, reduction = "pca", dims = 1:number_pcs,do.fast = T, k.seed = 10, check_duplicates = FALSE,perplexity = 30)
      
      
    }
    if (Calculate_MCPcounter) {
      #compute mcp scores
      if (normalization_method=="SCTransform") {
      mcp_scores <- MCPcounter.estimate(seurat_obj[["SCT"]]$data,featuresType = "HUGO_symbols")
      }
      if (normalization_method!="SCTransform") {
        mcp_scores <- MCPcounter.estimate(seurat_obj[["Spatial"]]$data,featuresType = "HUGO_symbols")
      }
      rownames(mcp_scores) <- make.names(rownames(mcp_scores))
      cell_types <- rownames(mcp_scores)
      for(cell_type in cell_types){
        seurat_obj <- AddMetaData(object= seurat_obj,metadata = mcp_scores[cell_type,],col.name = cell_type)
      }
      
    }
    if (generate_plots) {
      ####Plot The original lame + "nCount_Spatial", "nFeature_Spatial" and the spots removed...
      ##Get the sample name
      slide<-sub("\\/outs.*", "", raw_data_directory)
      slide<-sub(".*\\/", "", slide)
      
      he <- SpatialPlot(seurat_obj,repel = F,label = F,image.alpha=1,alpha = c(0,0)) + NoLegend() + ggtitle(slide)
      to_plot<-c("nCount_Spatial","nFeature_Spatial")
      removed=list(removed=setdiff(rownames(seurat_obj_before@meta.data),rownames(seurat_obj@meta.data)))
      spots_removed<-SpatialDimPlot(seurat_obj_before, cells.highlight =removed,facet.highlight = F, ncol = 1,pt.size.factor = 1.4)+ ggtitle(paste0("Spots removed: ", ncol(seurat_obj_before) - ncol(seurat_obj)))
      
      pl<- suppressMessages(lapply(to_plot,function(sign){
        SpatialPlot(seurat_obj,image.alpha = 0, pt.size.factor = 1.5,features = sign)+
          scale_fill_viridis_c(option='turbo') +
          DarkTheme() +
          hide_axis +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          theme(text=element_text(size=10))+
          theme(text=element_text(face = "bold"))+
          theme(legend.text=element_text(size=8))+theme(legend.text = element_text(angle = 90,margin = margin(8, 0, 0, 0, "pt")))
      }))
      pl <- append(list(he,spots_removed),pl)
      #save pdf
      ##Get the sample name
      slide<-sub("\\/outs.*", "", raw_data_directory)
      slide<-sub(".*\\/", "", slide)
      pdf(file = paste0(plots_path,"/",slide,"_",Sys.Date(),"_","Quality_spots_features.pdf"), width = 20, height = 6)
      print(patchwork::wrap_plots(pl,ncol = 4))  
      dev.off()  
      ###Plot MCP counter populations
      if (Calculate_MCPcounter) {
        ##Get the sample name
        slide<-sub("\\/outs.*", "", raw_data_directory)
        slide<-sub(".*\\/", "", slide)
        to_plot <- make.names(c(cell_types))
        
        he <- SpatialPlot(seurat_obj,repel = F,label = F,image.alpha=1,alpha = c(0,0)) + NoLegend() + ggtitle(slide)
        unsup_clustering <- suppressMessages(SpatialDimPlot(seurat_obj,label = T,label.size = 8)+NoLegend())
        pl<- suppressMessages(lapply(to_plot,function(sign){
          SpatialPlot(seurat_obj,image.alpha = 0, pt.size.factor = 1.5,features = sign)+
            scale_fill_viridis_c(option='turbo') +
            DarkTheme() +
            hide_axis +
            theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
            theme(text=element_text(size=10))+
            theme(text=element_text(face = "bold"))+
            theme(legend.text=element_text(size=8))
          #+theme(legend.text = element_text(angle = 90,margin = margin(8, 0, 0, 0, "pt")))
        }))
        pl <- append(list(he,unsup_clustering),pl)
        
        pdf(file = paste0(plots_path,"/",slide,"_",Sys.Date(),"_","MCPcounter.pdf"), width = 18, height = 8)
        print(patchwork::wrap_plots(pl,ncol = 6))  
        dev.off()  
        
      }
      
    }
    if (run_DE) {
      ##Find DE markers across clusters
      seurat_obj_markers <- FindAllMarkers(seurat_obj, min.pct = 0.25, logfc.threshold = 0.25,test.use = marker.test)
      return(list(seurat_obj,seurat_obj_markers))
    }else{
      return(seurat_obj)
    }
    
    
  }

library(SCP)
#library(SpaCET)
library(MCPcounter)
#####Importing starfysh deconvolution values and annotating each slide ############
##Import metaData
library(readr)
Starfysh_integrated_obs_all <- read_csv("Starfysh_integrated_obs_all.csv")##Path were you stored results

Starfysh_integrated_obs_all$barcode<-Starfysh_integrated_obs_all$...1
rownames(Starfysh_integrated_obs_all)<-Starfysh_integrated_obs_all$...1
Starfysh_integrated_obs_all<-Starfysh_integrated_obs_all[,-1]
Starfysh_integrated_obs_all<-as.data.frame(Starfysh_integrated_obs_all)

#########Charging, preprocessing and calculating GABA_synthesis + GABA_receptor for each slide
slides_id<-paste0("P_",seq(01,19,by=1)) ##This is the patient id as downloaded from GEO

meta_all<-data.frame()
for (i in 1:19) {
 ##Create a folder if it does not exist
  tmp_name<-slides_id[i]
  print(tmp_name)
   file.tmp<-paste0("MyoutsPath","/",tmp_name)
  if(!file.exists(file.tmp)){
    dir.create(file.path(file.tmp))
  }
  #setwd(file.tmp)  
setwd(file.tmp)
     visiumPath<-"Path_to_outsFolder" #Path to visium results
 Seurat_obj<-process_Visium_seurat(visiumPath,#Path to outs                                
                                    percent.mito_manual=30,
                                    min.genes=5,
                                     normalization_method="LogNormalize", #Options: SCTransform or NormalizeData
                                  nFeature=nFeature,
                                number_pcs=50,
                              use.jack=T,
                            resolution_clus=0.5,
                          Calculate_MCPcounter=T,
                        generate_plots=F,
                      plots_path="",
                    marker.test="roc",
                  only.pos = TRUE, 
                run_DE=F
 )
 
 ##Normalize using SCTransform but set variables to the same number so we can use the AddModuleScore without lossing genes due to variableFeatures
Seurat_obj <- SCTransform(Seurat_obj, verbose = T, 
                      #variable.features.n = nfeatures,
                      assay = "Spatial",vars.to.regress="percent.mito")

###Add name of sample
Seurat_obj$Sample<-tmp_name
     

  
##Calculate GABA signatures
intersect_sign <- intersect(Manual_signatures$GABA_synthesis,rownames(Seurat_obj))
print(length(intersect_sign))
#

Seurat_obj <- AddMetaData(Seurat_obj,apply(as.matrix(Seurat_obj[["SCT"]]$data[intersect_sign,]),2,mean),col.name = "GABA_synthesis_mean")
 
##Now for the receptor
 intersect_sign <- intersect(Manual_signatures$GABA_receptor,rownames(Seurat_obj))
print(length(intersect_sign))
#

Seurat_obj <- AddMetaData(Seurat_obj,apply(as.matrix(Seurat_obj[["SCT"]]$data[intersect_sign,]),2,mean),col.name = "GABA_receptor_mean")
  
  
####Add starfysh proportions from All_slides_integrated@meta.data
###Add other population proportions...


Seurat_obj@meta.data$barcode_starfysh<-paste0(gsub("[_].*","",rownames(Seurat_obj@meta.data)),"-",Seurat_obj$Sample)
 
table(Seurat_obj@meta.data$barcode_starfysh%in%Starfysh_integrated_obs_all$barcode)
#

#
sigs_starfysh<-colnames(Starfysh_integrated_obs_all)
sigs_starfysh<-sigs_starfysh[-1]#Take out barcode

for (j in c(sigs_starfysh)) {
  
  Seurat_obj@meta.data[,j]<-Starfysh_integrated_obs_all[match(Seurat_obj$barcode_starfysh,Starfysh_integrated_obs_all$barcode),j]
    
}
 




 ###Save object...
##Save Seurat object
saveRDS(Seurat_obj,file = paste0("Seurat_obj_",tmp_name,".rds"))


###Add metaData only to another object
meta_tmp<-Seurat_obj@meta.data
meta_all<-rbind(meta_all,meta_tmp)
rm(Seurat_obj)
gc()  
}

#Save meta_all
saveRDS(meta_all,file = "meta_all.rds")


########Calculations for panel h
nnls1=c("ccRCC_PT_like_starfysh","ccRCC_Proflammatory_like_starfysh","ccRCC_p53_EMT_starfysh","ccRCC_PI3K_starfysh",sigs_starfysh[1:20])
ma=c("GABA_synthesis_mean","GABA_receptor_mean")


##Correlation signatures...
meta<-meta_all

val_all = Reduce(rbind,lapply(ma, function(m){ #ma for modules
    Reduce(rbind,lapply(c(nnls1), function(ct){ ##Populations
        res = sapply(unique(meta$Sample), function(o){
            val = NA
            tryCatch(expr = {
                ct_nei = paste0(ct) ##Recover distance
                sub = meta[meta$Sample == o ,] ##subset to sample 
                model = lm(sub[,m] ~ sub[,ct_nei]) #module ~ pop_dist
                est = summary(model)$coefficients[2,1]
                pval = -log10(summary(model)$coefficients[2,4]) #pvalue
                val = pval*sign(est)
                return(val)
            }, error = function(e){return(NA)})
        })
        #res[abs(res)<pval_thresh] = 0
        #res = sign(res)
        #med = sum(res, na.rm = TRUE)
        #amp = sum(res == sign(med), na.rm = TRUE)
        med = median(res, na.rm = TRUE)
        amp = mean(sign(res) == sign(med), na.rm = TRUE)
        return(data.frame('module' = m, 'ct' = ct, 'med' = med, 'amp' = amp))
        return(data.frame('module' = m, 'ct' = ct, 'med' = med, 'amp' = amp))
    }))
}))

head(val_all)

###Med is the median of the (+-log10(pvalue)); amp(fraction samples)
##Negavite values for med means the correlation is negative;

#val_R$module = factor(val_R$module, levels = rev(ma[c(5:8,18,19,15,10,17,16,1:3,9,4,11:14,20)])) ##ordenate according to desire value...

val_all$ct = factor(val_all$ct, levels = rev(nnls1)) 
val_all$module=factor(val_all$module,levels = c("GABA_synthesis_mean","GABA_receptor_mean"))
val_all2<-val_all

val_all$highlight = (val_all$amp >= 0.5) & (abs(val_all$med) > 1.3)#Equivalent for pval<0.05
val_all$med[val_all$med > 4] = 4
val_all$med[val_all$med < -4] = -4
p2R = ggplot(val_all, aes(x=ct,y=module)) +
    geom_point(shape=21,aes(size=amp,fill=med,colour=highlight)) +
    scale_fill_gradient2(low = brewer_pal(palette = 'RdBu', direction = -1)(n = 3)[1],
                         mid = brewer_pal(palette = 'RdBu', direction = -1)(n = 3)[2], 
                         high = brewer_pal(palette = 'RdBu', direction = -1)(n = 3)[3]) +
    scale_colour_manual(breaks=c('FALSE','TRUE'),values=c('white','black')) + 
    theme_bw() +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          axis.text = element_text(size=14, colour = "black"),
          axis.text.y = element_text(size=12, colour = "black"),
          axis.text.x = element_text(size=12, colour = "black", angle = 90, hjust = 1),
          axis.title=element_blank(),
          panel.border = element_rect(size = 0.7, linetype = "solid", colour = "black"))
print(p2R+coord_flip())
p2R<-p2R+coord_flip()
###############Plot boxplot distribution of pvalues ##########
#For synthesis
  m<-ma[1]
  val = Reduce(rbind,lapply(nnls1, function(ct){
    res = sapply(unique(meta$Sample), function(o){
      val = NA
      tryCatch(expr = {
       # ct_dist = paste0(ct, '_dist')
        sub = meta[meta$Sample == o,]
        model = lm(sub[,m] ~ sub[,ct])
        est = summary(model)$coefficients[2,1]
        pval = -log10(summary(model)$coefficients[2,4])
        pval[pval > 10] = 10
        val = pval*sign(est)
        return(val)
      }, error = function(e){return(NA)})
    })
    return(data.frame('Sample' = unique(meta$Sample), 'ct' = ct, 'res' = res))
  }))
  val$ct = factor(val$ct, levels = rev(nnls1))
  val$Sample = factor(val$Sample)
  #val$RECIST_R_nR<-meta$RECIST_R_nR[match(val$Sample,meta$Sample)]
  p = ggplot(val, aes(x=ct,y=res)) +
    geom_boxplot(coef = 0, outlier.shape = NA) +
    #geom_violin() + 
    geom_jitter(colour = "gray") +
    geom_hline(yintercept = log10(0.05), linetype = 2) + 
    geom_hline(yintercept = -log10(0.05), linetype = 2) + 
    #scale_colour_manual(values = colors, breaks = names(colors)) +
    ggtitle(m) + 
    theme_bw() +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          axis.text = element_text(size=14, colour = "black"),
          axis.text.y = element_text(size=12, colour = "black"),
          axis.text.x = element_text(size=12, colour = "black", angle = 90, hjust = 1),
          axis.title=element_blank()) + NoLegend()
  print(p+coord_flip())
p_synthesis<-p+coord_flip()
#For receptor
 m<-ma[2]
  val = Reduce(rbind,lapply(nnls1, function(ct){
    res = sapply(unique(meta$Sample), function(o){
      val = NA
      tryCatch(expr = {
       # ct_dist = paste0(ct, '_dist')
        sub = meta[meta$Sample == o,]
        model = lm(sub[,m] ~ sub[,ct])
        est = summary(model)$coefficients[2,1]
        pval = -log10(summary(model)$coefficients[2,4])
        pval[pval > 10] = 10
        val = pval*sign(est)
        return(val)
      }, error = function(e){return(NA)})
    })
    return(data.frame('Sample' = unique(meta$Sample), 'ct' = ct, 'res' = res))
  }))
  val$ct = factor(val$ct, levels = rev(nnls1))
  val$Sample = factor(val$Sample)
  #val$RECIST_R_nR<-meta$RECIST_R_nR[match(val$Sample,meta$Sample)]
  p = ggplot(val, aes(x=ct,y=res)) +
    geom_boxplot(coef = 0, outlier.shape = NA) +
    #geom_violin() + 
    geom_jitter(colour = "gray") +
    geom_hline(yintercept = log10(0.05), linetype = 2) + 
    geom_hline(yintercept = -log10(0.05), linetype = 2) + 
    #scale_colour_manual(values = colors, breaks = names(colors)) +
    ggtitle(m) + 
    theme_bw() +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          axis.text = element_text(size=14, colour = "black"),
          axis.text.y = element_text(size=12, colour = "black"),
          axis.text.x = element_text(size=12, colour = "black", angle = 90, hjust = 1),
          axis.title=element_blank()) + NoLegend()
  print(p+coord_flip())
p_receptor<-p+coord_flip()+scale_x_discrete(position = "top")
#Plot together...

######################
##Plotting panel h
ggsave(ggarrange(p_synthesis+ggtitle(NULL),p2R+theme(axis.text.x =element_blank()),p_receptor+ggtitle(NULL),nrow = 1,align = "hv"),filename = "Figure3h.pdf",width = 14,height = 8)

Panel i

Code

##You need to run panel h to reproduce this figure
###This is an example for P_6
i=6
library(Seurat)
library(SPATA2)
spatial_correlation2<-function(spat_obj,var,threshold=0.5,cortest="pearson",colors_pal=c("firebrick1","gray72","palegreen","dodgerblue")){
    values <- NULL
    for(v in var){
      if(v %in% rownames(spat_obj)){
        spat_obj <- AddMetaData(spat_obj,
                                ifelse(spat_obj@assays$SCT@data[v,] > quantile(spat_obj@assays$SCT@data[v,],threshold),"high","low"),
                                col.name = paste0(v,"_cat"))
        values[[v]] <- spat_obj@assays$SCT@data[v,]
        
      }else{
        spat_obj <- AddMetaData(spat_obj,
                                ifelse(spat_obj@meta.data[,v] > quantile(spat_obj@meta.data[,v],threshold),"high","low"),
                                col.name = paste0(v,"_cat"))
        values[[v]] <- spat_obj@meta.data[,v]
      }
      
    }
    var_cat <- paste0(var,sep="_cat")
    spat_obj  <- AddMetaData(spat_obj,
                             factor(apply(spat_obj@meta.data[,var_cat],1,function(x) paste0(x,collapse = "_")),levels=rev(c("low_low",
                                                                                                                            "low_high",
                                                                                                                            "high_low",
                                                                                                                            "high_high"))),
                             col.name = paste0(var,collapse = "_"))
    get_cols <- which(table(spat_obj[[paste0(var,collapse = "_")]][,1])!=0)
    get_cols <- brewer.pal(4,"Spectral")[get_cols]
    # p1 <- SpatialDimPlot(spat_obj,
    #                      paste0(var,collapse = "_"),
    #                      cols = get_cols)
    p1 <- SpatialPlot(spat_obj,
                      paste0(var,collapse = "_"),
                      image.alpha = 0,
                      cols = get_cols)
    
    corres <- cor.test(values[[var[1]]],
                       values[[var[2]]],method = cortest,na.action="na.omit")
    
    
    table_cont <-spat_obj@meta.data[,var_cat]
    order_factor <- rev(c("high_high","high_low","low_high","low_low") )
    
    p2 <- plot_contingency_one_stack(data = table_cont,
                                     variables = colnames(table_cont),
                                     create_labels = F,
                                     order_factor = order_factor,
                                     palette = rev(brewer.pal(4,"Spectral")))
    
    p_cor <- ifelse(round(corres$p.value,digits = 3) == 0,"p < 10e-16",paste0("p = ",round(corres$p.value,digits = 3)))
    p_chisq <- ifelse(round(p2[[2]]$p.value,digits = 3) == 0,"p < 10e-16",paste0("p = ",round(p2[[2]]$p.value,digits = 3)))
    
    #p1 <- p1 + ggtitle(paste0(var,collapse = " "),
    #                   paste0("Correlation = ",round(corres$estimate,digits = 3)," ", p_cor))
    p1 <- p1 + ggtitle(paste0(var,collapse = " "),
                       paste0("Chisq.test: ",p_chisq, "\n Rho: ", round(corres$estimate,digits = 3)))+ scale_fill_manual(values=colors_pal) + DarkTheme() + hide_axis
    #####Change title to red and bold if significant chisq.test (<0.01)
    if (p2[[2]]$p.value<=0.01&corres$estimate>0) {
      p1 <- p1 +theme(plot.subtitle = element_text(color="red", face = "bold"))
    }
    
    #compute correlation plots
    v1 <- setNames(values[[1]],colnames(spat_obj))
    v2 <- setNames(values[[2]],colnames(spat_obj))
    
    corr <- plot_correlation(vect_x =v1,vect_y = v2,labls = var,method = cortest)
    p3 <- corr[[2]]+theme_classic()
    return(list("spatial"=p1,
                "cont_table"=p2,
                "spatial_cor"=p3))
  }

 Seurat_obj<-readRDS(file = paste0("Seurat_obj_",tmp_name,".rds"))
  
  
  
 #For number (i)
 he_p<- SpatialPlot(Seurat_obj,repel = F,label = F,image.alpha=1,alpha = c(0,0)) + NoLegend() + ggtitle(tmp_name)
  
  ####Signatures
  ###Extra for SPATA2
tmp_spata<-SPATA2::asSPATA2(
    object = Seurat_obj,
    sample_name = tmp_name,
    image_name = "slice1",
    spatial_method = "Visium"
)

##For number (iii, iv, v, vi, vii, viii, ix)
#Return
pl2<- suppressMessages(lapply(c("GABA_synthesis_mean","ccRCC_PT_like_starfysh","ccRCC_Proflammatory_like_starfysh","ccRCC_p53_EMT_starfysh","ccRCC_PI3K_starfysh","GABA_receptor_mean","Starfysh_Macrophage M2"),function(sign){
    SPATA2::plotSurface(object = tmp_spata,
            color_by = c(sign), display_image=T,#alpha_by = "SCGB3A1",
            smooth = T, normalize =F,smooth_span = 0.6,
            pt_size = 1.2) + 
    ggplot2::labs(title = sign)
  }))

#####for numbers (x, xi)

##Correlations...
var=list(c("GABA_synthesis_mean","ccRCC_PT_like_starfysh"),c("GABA_receptor_mean","Starfysh_Macrophage M2"))

#
#Grab Legend...to add at the end of each graph
tmp<-spatial_correlation2(spat_obj = Seurat_obj,var = c("GABA_synthesis_mean","ccRCC_PT_like_starfysh"),cortest  = "spearman")
tmp$spatial 
##Extract legend
library(ggpubr)
legend_spatial<-get_legend(tmp$spatial+ theme(legend.title = element_blank()) )
tmp


all_plots_cors<- lapply(var,function(v){
  tryCatch({
  res <- spatial_correlation2(spat_obj = Seurat_obj,
                             var = v,
                             cortest = "spearman")
  
  res$spatial <- res$spatial + hide_axis
  if(identical(v,var[length(var)][[1]])){
    res$spatial+ NoLegend()
  }else{
    res$spatial+ NoLegend()
  }
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
})

names(all_plots_cors)<-lapply(var, function(x){paste0(x[1],x[2])})

all_plots_cors